load c:\manoj\projects\plasma\HUA_PIU1 HUA_PIU_fday HUAPIU magnetic_frac;
data11 = load('c:\manoj\projects\plasma\STRICT_CRT_PLP_EOP_JT_09-17LT_268-298LON_KP.txt'); %15dec6 Claudia's new manual selecteion see email 
load c:\manoj\projects\plasma\Julia_PLP.mat;
fday_list_julia = floor(fday(:,1));
fday_julia_frac = fday(1,:) - floor(fday(1,:));
% 8.3.2007. Here the idea is to just plot 
PLP = PLP_new;% PLP_new is IRI corrected PLP data at 400 km (16Jan2007)
%Comment the previous line to use the old IRI-uncorrected PLP data
fday_list_julia = floor(fday(:,1));
bin_radius = 30/(60*24); %Radius of Julia data to be averaged in minutes
kp_limit = 3; 
delay_index = 1;
delay_sel = [13,14,15,16,17,18,19,20,21,22,23,24];
Kp_delay = delay_sel(delay_index);%Valid columns [10,13-24];
kp_delay_text = ['00';'01';'02';'03';'04';'05';'06';'07';'08';'09';'22';'MN'];
min_data = 10; %Minimum number of data pairs to be considered for correlation
robft = 1;
LT = [12 17]; %CHAMP LT limit 
LAG = [120]/(60*24); 
    
L = ~isnan(PLP(:,27)) & ~isnan(PLP(:,18))  & data11(:,Kp_delay) <=kp_limit & PLP(:,7) >= LT(1) & PLP(:,7) <= LT(2)...
          & data11(:,12) == 1;% & PLP(:,20) >= 1.5;%& PLP(:,29) >= 1.5;% & PLP(:,8) > -85 & PLP(:,8) < -70;
         IndData = PLP(L,:);
         counter = 0;
         
         for i = 1:size(IndData,1),
         fday_champ = tmjd(IndData(i,1),IndData(i,2),IndData(i,3),IndData(i,4),IndData(i,5),IndData(i,6),1);   
         iday = find(HUA_PIU_fday==floor(fday_champ));
         if ~isempty(iday),
            fday_champ_frac =  fday_champ - floor(fday_champ) - LAG(1);
            LLL = isnan(HUAPIU(iday,:));kkk = find(LLL == 0);
            if length(kkk) > 3,
             if fday_champ_frac > magnetic_frac(kkk(1)) & fday_champ_frac < magnetic_frac(kkk(end)),
                 counter=counter+1;
                  fday_champ_HUA(counter) = fday_champ;
                [trash left_bin] =  min(abs(magnetic_frac(kkk(1):kkk(end)) - (fday_champ_frac-bin_radius)));
                [trash right_bin] = min(abs(magnetic_frac(kkk(1):kkk(end)) - (fday_champ_frac+bin_radius)));
                left_bin  = kkk(1) + left_bin - 1;
                right_bin = kkk(1) + right_bin -1;
                HUA_data(counter) = nanmean(HUAPIU(iday,left_bin:right_bin));%<------- NORMAL
                %HUA_data(counter) = nanmean(HUAPIU(iday,860:1060)); % Daily
              
                Plp_ed(counter) = (IndData(i,18) + IndData(i,27))/2;
                Plp_lv(counter) = (IndData(i,17) + IndData(i,26))/2;
                Plp_ratio(counter) = (IndData(i,20) + IndData(i,29))/2;
                Plp_pt(counter) = Plp_ed(counter) -IndData(i,36);
                Plp_altitude(counter) = IndData(i,9);
             end;
            end;
         end;
         end;
        if counter >= min_data,
         LP = isnan(HUA_data);
        HUA_data(LP) = [];Plp_ed(LP) = [];Plp_lv(LP) = [];Plp_ratio(LP) = [];
        Plp_pt(LP) = [];fday_champ_HUA(LP) = [];
        counter = counter - sum(LP);
        end;
        if robft == 1 & counter >= min_data,
                temp_1 = rob_regress1(HUA_data,Plp_pt,Plp_ratio,Plp_lv,Plp_ed);
                HUA_data(temp_1) = [];
                Plp_ed(temp_1) = [];
                Plp_lv(temp_1) = [];
                Plp_ratio(temp_1) = [];
                Plp_pt(temp_1) = [];
                fday_champ_HUA(temp_1) = [];
                counter = length(HUA_data);
         end;
         lt_count = counter;
             [cof_ed,err_ed] = corrcoef(HUA_data,Plp_ed); 
             [cof_lv,err_lv] = corrcoef(HUA_data,Plp_lv); 
             [cof_rt,err_rt] = corrcoef(HUA_data,Plp_ratio); 
             [cof_pt,err_pt] = corrcoef(HUA_data,Plp_pt); 
             ED_matrix = cof_ed(2,1);
             LV_matrix = cof_lv(2,1);
             ED_error = err_ed(2,1);
             RT_matrix = cof_rt(2,1);
             RT_error = err_rt(2,1);
             PT_matrix = cof_pt(2,1);
             PT_error = err_pt(2,1);
             LV_error = err_lv(2,1);
         figure(1);    
         subplot(2,2,1);
         plot(HUA_data,Plp_ratio,'r.');
         title(['HUA ' sprintf(' Vs. CTR CC %4.2f',RT_matrix)]);
         subplot(2,2,2);
         plot(HUA_data,Plp_lv,'r.');
         title(['HUA ' sprintf(' Vs. L-value CC %4.2f',LV_matrix)]);
         subplot(2,2,3);
         plot(HUA_data,Plp_pt,'r.');
         title(['HUA ' sprintf(' Vs. PT CC%4.2f',PT_matrix)]);
         subplot(2,2,4);
         plot(HUA_data,Plp_ed,'r.');
         title(['HUA ' sprintf(' Vs. ED CC %4.2f',ED_matrix)]);
         
load c:\manoj\projects\plasma\plpf107hua HUA_data Plp_ratio Plp_lv Plp_ed Plp_pt fday_champ_HUA;
[a,b,c] = unique(floor(fday_champ_HUA));
load C:\Manoj\projects\plasma\solar_flux.mat
[a1,b1,c1] = intersect(a,fday_euvac);  
euvac_HUA = euvac(c1)';

corrcoef([HUA_data(b);Plp_ratio(b);Plp_lv(b);Plp_pt(b);Plp_ed(b);euvac_HUA]')
        figure(2);
         subplot(2,2,1);
         plot(euvac_HUA,Plp_ratio(b),'r.');
         c = corrcoef(euvac_HUA,Plp_ratio(b));
         title(sprintf(' P Vs. CTR CC %4.2f',c(2,1)));
         subplot(2,2,2);
         plot(euvac_HUA,Plp_lv(b),'r.');
         c = corrcoef(euvac_HUA,Plp_lv(b));
         title(sprintf(' P Vs. L-value CC %4.2f',c(2,1)));
         subplot(2,2,3);
         plot(euvac_HUA,Plp_pt(b),'r.');
         c = corrcoef(euvac_HUA,Plp_pt(b));
         title(sprintf(' P Vs. PT CC%4.2f',c(2,1)));
         subplot(2,2,4);
         plot(euvac_HUA,Plp_ed(b),'r.');
         c = corrcoef(euvac_HUA,Plp_ed(b));
         title(sprintf(' P Vs. ED CC %4.2f',c(2,1)));


         
load c:\manoj\projects\plasma\plpf107hua HUA_data Plp_ratio Plp_lv Plp_ed Plp_pt fday_champ_HUA;
[a,b,c] = unique(floor(fday_champ_HUA));
load C:\Manoj\projects\plasma\solar_flux.mat
[a1,b1,c1] = intersect(a,fday_euvac);  
euvac_HUA = euvac(c1)';
         
ctr_euvac = Plp_ratio(b);
pt_euvac = Plp_pt(b);
HUA_euvac = HUA_data(b);
lv_euvac = Plp_lv(b);
ed_euvac = Plp_ed(b);

pp=1;
%f107bin = [70:30:200];
 f107bin = [70:30:100];
for i = 1:length(f107bin)-1,
    L = euvac_HUA >= f107bin(i) & euvac_HUA < f107bin(i+1);
    subplot(2,2,1);
    plot(HUA_euvac(L),ctr_euvac(L),'k.','MarkerSize',16);hold on;
    c=corrcoef(HUA_euvac(L),ctr_euvac(L));
    [c1,sfit] = polyfit(HUA_euvac(L),ctr_euvac(L),1);
    m = polyval(c1,[min(HUA_euvac(L)),max(HUA_euvac(L))]);
    plot([min(HUA_euvac(L)),max(HUA_euvac(L))],m,'r-','LineWidth',2);hold off;
    axis([-50 150 1 4]);
    text(100,1.5,sprintf('%5.2e %5.2f',c1(1),c1(2)));fprintf('%f %f\n',c1(1),c1(2));
    xlabel('HUA_{\DeltaH} - PIU_{\DeltaH}');
    ylabel('Ctrest-Trough Ratio');
    title(sprintf('EUVAC from %d to %d, CC %5.2f',f107bin(i),f107bin(i+1),c(2,1)));
subplot(2,2,2);
    plot(HUA_euvac(L),lv_euvac(L),'k.','MarkerSize',16);hold on;
    [c1,sfit] = polyfit(HUA_euvac(L),lv_euvac(L),1);
    m = polyval(c1,[min(HUA_euvac(L)),max(HUA_euvac(L))]);
    plot([min(HUA_euvac(L)),max(HUA_euvac(L))],m,'r-','LineWidth',2);hold off;
    axis([-50 150 1 1.3]);
    text(100,1.05,sprintf('%5.2e %5.2f',c1(1),c1(2)));fprintf('%f %f\n',c1(1),c1(2));
    c=corrcoef(HUA_euvac(L),lv_euvac(L));
    title(sprintf('EUVAC from %d to %d, CC %5.2f',f107bin(i),f107bin(i+1),c(2,1)));
    xlabel('HUA_{\DeltaH} - PIU_{\DeltaH}');
    ylabel('L-Value');


subplot(2,2,3);
    plot(HUA_euvac(L),pt_euvac(L),'k.','MarkerSize',16);hold on;
    [c1,sfit] = polyfit(HUA_euvac(L),pt_euvac(L),1);
    m = polyval(c1,[min(HUA_euvac(L)),max(HUA_euvac(L))]);
    plot([min(HUA_euvac(L)),max(HUA_euvac(L))],m,'r-','LineWidth',2);hold off;
    axis([-50 150 0 2e6]);
    text(80,1e6,sprintf('%5.2e %5.2e',c1(1),c1(2)));fprintf('%f %f\n',c1(1),c1(2));

        c=corrcoef(HUA_euvac(L),pt_euvac(L));
    title(sprintf('EUVAC from %d to %d, CC %5.2f',f107bin(i),f107bin(i+1),c(2,1)));
    xlabel('HUA_{\DeltaH} - PIU_{\DeltaH}');
    ylabel('Peak-Trough');


subplot(2,2,4);
    plot(HUA_euvac(L),ed_euvac(L),'k.','MarkerSize',16);  
    c=corrcoef(HUA_euvac(L),ed_euvac(L));hold on;
    [c1,sfit] = polyfit(HUA_euvac(L),ed_euvac(L),1);
    m = polyval(c1,[min(HUA_euvac(L)),max(HUA_euvac(L))]);
    plot([min(HUA_euvac(L)),max(HUA_euvac(L))],m,'r-','LineWidth',2);hold off;
    axis([-50 150 0 4e6]);
    text(80,1e6,sprintf('%5.2e %5.2e',c1(1),c1(2)));fprintf('%f %f\n',c1(1),c1(2));
        
    title(sprintf('EUVAC from %d to %d, CC %5.2f',f107bin(i),f107bin(i+1),c(2,1)));
    xlabel('HUA_{\DeltaH} - PIU_{\DeltaH}');
    ylabel('mean peak electron density');

    pause;
end;
    
    
